psychencode_figures_cleaned

Packages etc.

library(dplyr)
library(tidyverse)
library(edgeR)
library(markerGeneProfile) 
library(matrixStats)
library(cowplot) 
library(broom)
library(knitr)
library(ggpubr)
library(ggrepel)
library(patchwork)
library(ggsignif)
library(modelr)
library(ggbeeswarm)
library(lemon)
theme_set(theme_classic2())
#Colour palette
cbPalette = c("#56B4E9", "#009E73","#E69F00", "#0072B2", "#D55E00", "#CC79A7","#000000","#F0E442")

Loading Data

#Cell type proportions for each data set merged with metadata - ready for modelling 
mgp_estimations = read_csv( "/external/rprshnas01/kcni/dkiss/cell_prop_psychiatry/data/psychencode_mgp_estimations.csv") 
acc_estimations = read_csv("/external/rprshnas01/kcni/dkiss/cell_prop_psychiatry/data/cmc_acc_mgp_estimations.csv") %>% filter(primaryDiagnosis %in% c("control", "Schizophrenia"))
#Factorize newStudy for aesthetic purposes (this orderl ooks good on graphs)
mgp_estimations = mgp_estimations %>% filter(ageDeath >= 15)
mgp_estimations$newStudy = mgp_estimations$newStudy %>% factor(levels = c("Pitt", "GVEX", "NIMH_HBCC", "LIBD_szControl", "MSSM", "Penn"))
mgp_estimations_long = mgp_estimations %>% pivot_longer(cols = Astrocyte:VLMC, names_to = 'cell_type', values_to = 'rel_prop') %>%
  mutate(cell_class = case_when(
    cell_type %in% c( "PVALB", "SST", "VIP", "LAMP5", "PAX6") ~ "Inhibitory",
    cell_type %in% c("IT", "L4.IT", "L5.ET", "L5.6.IT.Car3", "L5.6.NP", "L6.CT", "L6b") ~ "Excitatory",
    TRUE ~ "Non-Neuronal"
  ))
mgp_estimations_long$cell_type = mgp_estimations_long$cell_type %>% factor(levels = c("PVALB", "SST", "VIP", "LAMP5", "PAX6", "IT", "L4.IT", "L5.ET", "L5.6.IT.Car3", "L5.6.NP", "L6.CT", "L6b", "Astrocyte", "Endothelial", "Microglia", "Oligodendrocyte", "OPC", "Pericyte", "VLMC"))
mgp_estimations_long$cell_class = mgp_estimations_long$cell_class %>% factor(levels = c("Inhibitory", "Excitatory", "Non-Neuronal"))

acc_estimations_long = acc_estimations %>% pivot_longer(cols = Astrocyte:VLMC, names_to = 'cell_type', values_to = 'rel_prop') %>%
  mutate(cell_class = case_when(
    cell_type %in% c("PVALB", "SST", "VIP", "LAMP5", "PAX6") ~ "Inhibitory",
    cell_type %in% c("IT", "L4.IT", "L5.ET", "L5.6.IT.Car3", "L5.6.NP", "L6.CT", "L6b") ~ "Excitatory",
    TRUE ~ "Non-Neuronal"
  ))
acc_estimations_long$cell_type = acc_estimations_long$cell_type %>% factor(levels = c("PVALB", "SST", "VIP", "LAMP5", "PAX6", "IT", "L4.IT", "L5.ET", "L5.6.IT.Car3", "L5.6.NP", "L6.CT", "L6b", "Astrocyte", "Endothelial", "Microglia", "Oligodendrocyte", "OPC", "Pericyte", "VLMC"))
acc_estimations_long$cell_class = acc_estimations_long$cell_class %>% factor(levels = c("Inhibitory", "Excitatory", "Non-Neuronal"))

Load snRNAseq proportions

psychencode_snCTPs = read.csv('/external/rprshnas01/netdata_kcni/stlab/Xiaolin/cell_deconv_data/PsychEncode_label_transferred_snCTP.csv')
snCTP_metadata = read.csv("~/cell_prop_psychiatry/data/psychencode_label_transferred_pseudobulk_metadata.csv")[,-1] %>%
  group_by(ID) %>%
  mutate(total_cells = sum(num_cells)) %>%
  ungroup() 
#Add log10 number of cells variable
snCTP_metadata$log10_num_cells = log10(snCTP_metadata$num_cells)
#Add age group variable
snCTP_metadata =  snCTP_metadata %>%
  mutate(age_group = case_when(
    Age < 70 ~ "Under_70", Age >= 70 ~ "Over_70")) 
snCTP_metadata$age_group = snCTP_metadata$age_group %>% factor(levels = c("Under_70", "Over_70"))
psychencode_snCTPs = left_join(snCTP_metadata, psychencode_snCTPs) 

ling_snCTPs = read.csv("~/cell_prop_psychiatry/data/ling_snCTP_estimations")

Cohort demographic summary

summary_df = mgp_estimations %>%
  mutate(individualIDSource = factor(newStudy, levels = c("Pitt", "GVEX", "NIMH_HBCC", "LIBD_szControl", "MSSM", "Penn"))) %>%
    group_by(newStudy) %>%
    summarise(
        CohortName = first(newStudy),
        PercentMale = mean(reportedGender == "male", na.rm = TRUE) * 100,
        AvgAge = mean(ageDeath, na.rm = TRUE),
        AgeSD = sd(ageDeath, na.rm = TRUE),
        NumCases = sum(primaryDiagnosis == "Schizophrenia", na.rm = TRUE),
        NumControls = sum(primaryDiagnosis == "control", na.rm = TRUE)
    )

Residualize MGPs based on linear model

mgp_estimations_long = mgp_estimations_long %>%
  group_by(cell_type) %>%
  mutate(resid = resid(lm(scale(rel_prop) ~ scale(RIN) + scale(ageDeath) + scale(PMI) + 
                reportedGender))) %>%
  ungroup()

Model rCTPs based on diagnosis and covariates seperately in each cohort

combined_lms = mgp_estimations_long %>%
  group_by(newStudy, cell_type) %>%
  do(tidy(lm(scale(rel_prop) ~ scale(RIN) + scale(ageDeath) + scale(PMI) +
               reportedGender + primaryDiagnosis, data = .))) %>%
  ungroup() %>%
  mutate(padj = p.adjust(p.value, method = 'fdr'),
         term = recode(term,
                       `(Intercept)` = "Intercept",
                       `reportedGendermale` = "gender:Male",
                       `primaryDiagnosisSchizophrenia` = "SCZ",
                       `scale(ageDeath)` = "Age",
                       `scale(PMI)` = "PMI",
                       `scale(RIN)` = "RIN")) %>%
  merge(unique(mgp_estimations_long[c('dataset', 'newStudy')]), by.x = "newStudy") %>%
  merge(unique(mgp_estimations_long[c('cell_class', 'cell_type')]), by.x = "cell_type") %>%
  mutate(newStudy = factor(newStudy, levels = c("Pitt", "GVEX", "NIMH_HBCC", "LIBD_szControl", "MSSM", "Penn")))

Interaction Models: As outlined above, there appears to be an interaction between age x diagnosis with respect to SST cell proportion. In other words, SST cell proportion is generally lower as ageDeath increases, but declines faster with age in controls. To analyze this more clearly, we create new models (interaction_lms) to visualize this effect. These models include new_study as a covariate, and age x diagnosis as an interactor.

interaction_lms = mgp_estimations_long %>%
  group_by(newStudy, cell_type) %>%
  do(tidy(lm(scale(rel_prop) ~ scale(RIN) + scale(ageDeath) + scale(PMI) +
               reportedGender + primaryDiagnosis + scale(ageDeath) * primaryDiagnosis, data = .))) %>%
  ungroup() %>%
  mutate(padj = p.adjust(`p.value`, method = 'fdr'),
         term = recode(term,
                       `(Intercept)` = "Intercept",
                       `reportedGendermale` = "gender:Male",
                       `primaryDiagnosisSchizophrenia` = "SCZ",
                       `scale(ageDeath):primaryDiagnosisSchizophrenia` = "Age x Diagnosis",
                       `scale(ageDeath)` = "Age",
                       `scale(PMI)` = "PMI",
                       `scale(RIN)` = "RIN")) %>%
  merge(unique(mgp_estimations_long[c('dataset', 'newStudy')]), by.x = "newStudy") %>%
  merge(unique(mgp_estimations_long[c('cell_class', 'cell_type')]), by.x = "cell_type") %>%
  mutate(newStudy = factor(newStudy, levels = c("Pitt", "GVEX", "NIMH_HBCC", "LIBD_szControl", "MSSM", "Penn")))

Interaction mega-analysis w/ all data

mega_lms = mgp_estimations_long %>% 
  group_by(cell_type) %>%
  do(tidy(lm(scale(rel_prop) ~ scale(RIN) + scale(ageDeath) + scale(PMI) +
             reportedGender + primaryDiagnosis + newStudy + scale(ageDeath)*primaryDiagnosis, data = .))) %>%
  ungroup() %>% 
  mutate(padj = p.adjust(`p.value`, method = 'fdr'),
         term = recode(term,
                       `(Intercept)` = "Intercept",
                       `reportedGendermale` = "gender:Male",
                       `primaryDiagnosisSchizophrenia` = "SCZ",

                       `scale(ageDeath)` = "Age",
                       `scale(ageDeath):primaryDiagnosisSchizophrenia` = "Age x Diagnosis",
                       `scale(PMI)` = "PMI", 
                       `newStudyNIMH_HBCC` = "NIMH_HBCC",
                       `newStudyLIBD_szControl` = "LIBD_szControl",
                       `newStudyPitt` = "Pitt",
                       `newStudyMSSM` = "MSSM",
                       `newStudyPenn` = "Penn",
                       `scale(RIN)` = "RIN")) %>%
  merge(unique(mgp_estimations_long[c('cell_class', 'cell_type')]), by.x = "cell_type") 

Functions to generate age range bins

# Create function to map age to age range
age_to_range = function(age) {
    if(age >= 15 & age < 30) {
    return("15-29")
  } else {
    lower_bound = floor((age - 30) / 10) * 10 + 30
    upper_bound = lower_bound + 9
    return(paste0(lower_bound, "-", upper_bound))
  }
}
# Create function to bin ages to >=70 and <70
classify_ages = function(ages) {
  classifications = ifelse(ages <= 70, "Under 70", "Over 70")
  return(classifications)
}

Create age groups from age 20-90 in increments of 10 years. Plot meganalysis beta coeffs against age group.

mgp_estimations_long$age_range = sapply(mgp_estimations_long$ageDeath, age_to_range)

mega_lms_by_age = mgp_estimations_long %>% 
  group_by(age_range, cell_type) %>%
  do(tidy(lm(scale(rel_prop) ~ scale(RIN) + scale(PMI) +
             reportedGender + primaryDiagnosis + newStudy, data = .))) %>%
  ungroup() %>% 
  mutate(padj = p.adjust(`p.value`, method = 'fdr'),
         term = recode(term,
                       `(Intercept)` = "Intercept",
                       `reportedGendermale` = "gender:Male",
                       `primaryDiagnosisSchizophrenia` = "SCZ",
                       `scale(ageDeath)` = "Age",
                       `scale(PMI)` = "PMI", 
                       `newStudyNIMH_HBCC` = "NIMH_HBCC",
                       `newStudyPitt` = "Pitt",
                       `newStudyMSSM` = "MSSM",
                       `newStudyPenn` = "Penn",
                       `scale(RIN)` = "RIN")) %>%
  merge(unique(mgp_estimations_long[c('cell_class', 'cell_type')]), by.x = "cell_type") 

Bulk RNAseq Meganalysis based on age +/- 70

mgp_estimations_long$age_class <- sapply(mgp_estimations_long$ageDeath, classify_ages) %>% factor(levels = c("Under 70", "Over 70"))
mega_lms_by_age_class = mgp_estimations_long %>% 
  group_by(age_class, cell_type) %>%
  do(tidy(lm(scale(rel_prop) ~ scale(RIN) + scale(PMI) +
             reportedGender + newStudy + primaryDiagnosis, data = .))) %>%
  ungroup() %>% 
  mutate(padj = p.adjust(`p.value`, method = 'fdr'),
         term = recode(term,
                       `(Intercept)` = "Intercept",
                       `reportedGendermale` = "gender:Male",
                       `primaryDiagnosisSchizophrenia` = "SCZ",
                       `scale(ageDeath)` = "Age",
                       `scale(PMI)` = "PMI", 
                       `newStudyNIMH_HBCC` = "NIMH_HBCC",
                       `newStudyPitt` = "Pitt",
                       `newStudyMSSM` = "MSSM",
                       `newStudyPenn` = "Penn",
                       `scale(RIN)` = "RIN")) %>%
  merge(unique(mgp_estimations_long[c('cell_class', 'cell_type')]), by.x = "cell_type") 

#same but stratify by cohort
lms_by_age_class = mgp_estimations_long %>% 
  group_by(age_class, cell_type, newStudy) %>%
  do(tidy(lm(scale(rel_prop) ~ scale(RIN) + scale(PMI) +
             reportedGender + primaryDiagnosis, data = .))) %>%
  ungroup() %>% 
  mutate(padj = p.adjust(`p.value`, method = 'fdr'),
         term = recode(term,
                       `(Intercept)` = "Intercept",
                       `reportedGendermale` = "gender:Male",
                       `primaryDiagnosisSchizophrenia` = "SCZ",
                       `scale(ageDeath)` = "Age",
                       `scale(PMI)` = "PMI", 
                       `newStudyNIMH_HBCC` = "NIMH_HBCC",
                       `newStudyPitt` = "Pitt",
                       `newStudyMSSM` = "MSSM",
                       `newStudyPenn` = "Penn",
                       `scale(RIN)` = "RIN")) %>%
  merge(unique(mgp_estimations_long[c('cell_class', 'cell_type')]), by.x = "cell_type") 

Single Cell Meganalysis with interaction model

snCTPs_long = psychencode_snCTPs %>%
  pivot_longer(cols = Astrocyte:VIP, names_to = 'cell_type', values_to = 'rel_prop',)  %>%
  mutate(cell_class = case_when(
    cell_type %in% c("LAMP5", "PAX6", "PVALB", "SST", "VIP") ~ "Inhibitory",
    cell_type %in% c("IT", "L4.IT", "L5.ET", "L5.6.IT.Car3", "L5.6.NP", "L6.CT", "L6b") ~ "Excitatory",
    TRUE ~ "Non-Neuronal"
  )) %>%  
  distinct(cell_type, unique_donor_ID, .keep_all = T)


sn_mega_lms = snCTPs_long %>% filter(total_cells >= 500) %>%
  group_by(cell_type) %>%
  do(tidy(lm(scale(rel_prop) ~ scale(PMI) + Gender + Institution + Phenotype + scale(Age) * Phenotype, data = .))) %>%
  ungroup() %>%
  mutate(padj = p.adjust(`p.value`, method = 'fdr')) %>%
  # add cell class labels
  mutate(cell_class = case_when(
    cell_type %in% c("LAMP5", "PAX6", "PVALB", "SST", "VIP") ~ "Inhibitory",
    cell_type %in% c("IT", "L4.IT", "L5.ET", "L5.6.IT.Car3", "L5.6.NP", "L6.CT", "L6b") ~ "Excitatory",
    TRUE ~ "Non-Neuronal"
  )) %>%
  mutate(term = recode(term,
                       `(Intercept)` = "Intercept",
                       `reportedGendermale` = "gender:Male",
                       `PhenotypeSZ` = "SCZ",
                       `primaryDiagnosisBipolar Disorder` = "BP",
                       `DxMDD` = "MDD",
                       `scale(Age)` = "Age",
                       `scale(PMI)` = "PMI", 
                       `newStudyNIMH_HBCC` = "NIMH_HBCC",
                       `newStudyPitt` = "Pitt",
                       `newStudyMSSM` = "MSSM",
                       `newStudyPenn` = "Penn",
                       `scale(RIN)` = "RIN",
                       `PhenotypeSZ:scale(Age)` = "Interaction"
))

Single Cell Meganalysis based on age +/- 70

sn_mega_lms_by_age = snCTPs_long %>% filter(total_cells >= 500) %>%
  group_by(cell_type, age_group) %>%
  do(tidy(lm(scale(rel_prop) ~ scale(PMI) + Gender + Institution + Phenotype, data = .))) %>%
  ungroup() %>%
  mutate(padj = p.adjust(`p.value`, method = 'fdr')) %>%
  mutate(term = recode(term,
                       `(Intercept)` = "Intercept",
                       `reportedGendermale` = "gender:Male",
                       `PhenotypeSZ` = "SCZ",
                       `primaryDiagnosisBipolar Disorder` = "BP",
                       `DxMDD` = "MDD",
                       `scale(ageDeath)` = "Age",
                       `scale(PMI)` = "PMI", 
                       `newStudyNIMH_HBCC` = "NIMH_HBCC",
                       `newStudyPitt` = "Pitt",
                       `newStudyMSSM` = "MSSM",
                       `newStudyPenn` = "Penn",
                       `scale(RIN)` = "RIN"))

#again but stratified by cohort
sn_lms_by_age = snCTPs_long %>% filter(total_cells >= 500) %>%
  group_by(cell_type, age_group, Institution) %>%
  do(tidy(lm(scale(rel_prop) ~ scale(PMI) + Gender + Phenotype, data = .))) %>%
  ungroup() %>%
  mutate(padj = p.adjust(`p.value`, method = 'fdr')) %>%
  mutate(term = recode(term,
                       `(Intercept)` = "Intercept",
                       `reportedGendermale` = "gender:Male",
                       `PhenotypeSZ` = "SCZ",
                       `primaryDiagnosisBipolar Disorder` = "BP",
                       `DxMDD` = "MDD",
                       `scale(ageDeath)` = "Age",
                       `scale(PMI)` = "PMI", 
                       `newStudyNIMH_HBCC` = "NIMH_HBCC",
                       `newStudyPitt` = "Pitt",
                       `newStudyMSSM` = "MSSM",
                       `newStudyPenn` = "Penn",
                       `scale(RIN)` = "RIN"))

————////FIGURES////————

Figure 1: Illustrating the association between cell prop and diagnosis and cohort specific differences

significant_combined = subset(combined_lms, padj < 0.1 & cell_type %in% c("SST", "PVALB", "VIP") & term == "SCZ")
significant_combined$asterisks <- ifelse(significant_combined$padj <= 0.01, "***",
                                         ifelse(significant_combined$padj <= 0.05, "**", "*"))
significant_combined$vjust <- ifelse(significant_combined$estimate >= 0, -3, 7)

figure_1b = combined_lms %>% 
  filter(term == 'SCZ', cell_type %in% c("SST", "PVALB", "VIP")) %>% 
  mutate(cell_type = fct_reorder(cell_type, estimate)) %>% 
  ggplot(aes(x = newStudy, y = estimate, fill = newStudy)) + 
  geom_hline(yintercept = 0) + 
  geom_bar(stat = "identity", position = "dodge") + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error)) + 
  ylab('Beta Coefficient') + 
  xlab('') +
  scale_fill_manual(values = cbPalette) +
  geom_text(data = significant_combined,
            aes(newStudy, estimate, label = asterisks), 
            vjust = significant_combined$vjust) +
  theme(
    axis.text.x = element_blank(),
    axis.text.y = element_text(size = 15),
    axis.title.x = element_text(vjust = 6, size = 15),
    axis.title.y = element_text(size = 16),
    strip.text.x = element_text(size = 14),
    axis.line.x = element_blank(), 
    axis.ticks.x = element_blank()) +
  facet_grid(~cell_type, drop = T, scale = "free_x", space = "free") +
  guides(fill = guide_none()) +
  scale_y_continuous(limits = c(-0.9, 0.5))

# Figure 1c: Residualized MGPs vs. Diagnosis
figure_1c <- mgp_estimations_long %>%
  filter(cell_type %in% c("SST", "PVALB", "VIP")) %>%
  ggplot(aes(x = primaryDiagnosis, y = resid)) +
  geom_boxplot(outlier.shape = NA, aes(fill = newStudy, alpha = 0.5), show.legend = F, notch = T) +
  geom_beeswarm(size = 1, alpha = 0.3, aes(colour = newStudy), show.legend = F) +
  ylab('Residualized CTPs') + 
  xlab('Diagnosis') + 
  scale_x_discrete(labels = c("C", "SCZ")) +
  scale_colour_manual(values = cbPalette) +
  scale_fill_manual(values = cbPalette) +
  theme(
    axis.text.x = element_text(size = 15),
    axis.text.y = element_text(size = 15),
    axis.title.x = element_text(vjust = -2, size = 15),
    axis.title.y = element_text(size = 15),
    strip.text.x = element_text(size = 14),
    strip.text.y = element_text(size = 14)) +
  guides(fill = FALSE) +
  facet_grid(rows = vars(cell_type), cols = vars(newStudy), drop = T, scale = "free", switch = "y", space = "free") +
  geom_signif(comparisons = list(c("control", "Schizophrenia")), map_signif_level = TRUE, , test = wilcox.test) +
  scale_y_continuous(limits = c(-3.5, 3.5)) 

Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as of ggplot2 3.3.4.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

#Figure 1d - Hisotgram of cohort ages
figure_1d = mgp_estimations %>% 
    ggplot(aes(x = ageDeath)) +
    geom_histogram(aes(fill = newStudy)) +
    facet_wrap(~ newStudy, nrow = 6, scales = "free_y") +
    labs(x = "ageDeath", y = "Frequency") +
    scale_fill_manual(values = c("#56B4E9", "#009E73","#E69F00", "#0072B2", "#D55E00", "#CC79A7")) + 
  scale_y_continuous(n.breaks = 4) +
  guides(fill = guide_none()) + 
  theme(panel.grid = element_blank(),
        strip.background = element_blank(),
        strip.text = element_text(size = rel(0.9), face = "plain"),
    axis.text.x = element_text(size = 13,),
    axis.text.y = element_text(size = 13),
    axis.title.x = element_text( size = 15),
    axis.title.y = element_text(size = 16),
    strip.text.x = element_text(size = 14))

# Combine Figure 1
figure_1_top = (figure_1b | figure_1d) + plot_layout(widths = c(4,2))
figure_1_bottom = figure_1c
  
figure_1 =  figure_1_top / figure_1_bottom + plot_annotation(tag_levels = 'A') + plot_layout(heights = c(1, 1)) & theme(plot.tag = element_text(size = 25))
figure_1

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Warning: Removed 11 rows containing non-finite outside the scale range (`stat_boxplot()`).
Warning: Removed 11 rows containing non-finite outside the scale range (`stat_signif()`).
Warning: Removed 11 rows containing missing values or values outside the scale range (`geom_point()`).

#Save to plots folder
output_filename <- "~/cell_prop_psychiatry/plots/figure_1.jpg"

# Use ggsave to save the plot as a JPG image
ggsave(output_filename, figure_1, width = 11, height = 13, dpi = 500)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Warning: Removed 11 rows containing non-finite outside the scale range (`stat_boxplot()`).
Warning: Removed 11 rows containing non-finite outside the scale range (`stat_signif()`).
Warning: Removed 11 rows containing missing values or values outside the scale range (`geom_point()`).

Figure 2: Introducing and presenting an interaction between age and diagnosis w/rt cell proportion.

# Plot interaction beta values from meganalysis
significant_mega = subset(mega_lms, padj <= 0.1 & term == "Age x Diagnosis" & cell_type %in% c("Astrocyte", "IT", "L5.6.NP", "LAMP5", "Microglia", "Oligodendrocyte", "OPC", "Pericyte", "PVALB", "SST", "VIP" ))
significant_mega$asterisks = ifelse(significant_mega$padj <= 0.01, "***",
                                   ifelse(significant_mega$padj <= 0.05, "**", "*"))
significant_mega$vjust = ifelse(significant_mega$estimate >= 0, -3, 3)

# Figure 2a: Results from interaction meganalysis stratified by cell class - plot interaction beta coeff
mega_lms$cell_type = mega_lms$cell_type %>% factor(levels = c("PVALB", "SST", "VIP", "LAMP5", "PAX6", "IT", "L4.IT", "L5.ET", "L5.6.IT.Car3", "L5.6.NP", "L6.CT", "L6b", "Astrocyte", "Endothelial", "Microglia", "Oligodendrocyte", "OPC", "Pericyte", "VLMC"))
figure_2a = mega_lms %>%  filter(cell_type %in% c("Astrocyte", "IT", "L5.6.NP", "LAMP5", "Microglia", "Oligodendrocyte", "OPC","Pericyte", "PVALB", "SST", "VIP" )) %>%
  filter(term %in% 'Age x Diagnosis') %>% 
  ggplot(aes(x = cell_type, y = estimate)) + 
  geom_hline(yintercept = 0) + 
  geom_bar(stat = "identity", position = "dodge") + 
  facet_grid(~cell_class, drop = T, scale = "free_x", space = "free") +
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error)) + 
  ylab('Interaction Coefficient') + 
  xlab('Cell Type') + 
  geom_text(data = significant_mega,
            aes(cell_type, estimate, label = asterisks), 
            vjust = significant_mega$vjust) +
  theme(
    axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1, size = 12), 
    axis.text.y = element_text(size = 15),
    axis.title.x = element_text(vjust = 8, size = 15),
    axis.title.y = element_text(size = 15),
    strip.text.x = element_text(size = 14)) +
  guides(fill = FALSE) +
  scale_y_continuous(limits = c(-0.15, 0.32))

# Figure 2b: Plot beta values from meganalysis in 10-year age bins for SST, PVALB, VIP
significant_mega_age = subset(mega_lms_by_age, padj <= 0.1 & term == "SCZ" & cell_type %in% c("SST", "PVALB", "VIP"))
significant_mega_age$asterisks = ifelse(significant_mega_age$padj <= 0.01, "***",
                                         ifelse(significant_mega_age$padj <= 0.05, "**", "*"))
significant_mega_age$vjust = ifelse(significant_mega_age$estimate >= 0, -2, 3)

figure_2b = mega_lms_by_age %>% 
  filter(term %in% 'SCZ' & cell_type %in% c("SST", "PVALB", "VIP")) %>% 
  mutate(cell_type = fct_reorder(cell_type, estimate)) %>% 
  ggplot(aes(x = age_range, y = estimate)) + 
  geom_hline(yintercept = 0) + 
  geom_bar(stat = "identity", position = "dodge") + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error)) + 
  ylab('Beta Coefficient') + 
  xlab('\n \n Age Range') + 
  geom_text(data = significant_mega_age,
            aes(age_range, estimate, label = asterisks), 
            vjust = significant_mega_age$vjust) +
   theme(
    axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1, size = 13), 
    axis.text.y = element_text(size = 15),
    axis.title.x = element_text(vjust = 8, size = 15),
    axis.title.y = element_text(size = 15),
    strip.text.x = element_text(size = 14)) +
  facet_wrap(~cell_type, drop = T, scale = "free_x") +
  guides(fill = FALSE) +
  scale_y_continuous(limits = c(-1.0, 1.4))


# Figure 2c: Visualizing interaction between age and diagnosis by plotting rel_prop vs. age in each diagnosis 
figure_2c = mgp_estimations_long %>% filter(ageDeath < 90) %>%
  filter(cell_type %in% c("SST", "PVALB", "VIP")) %>%
  ggplot(aes(x = ageDeath, y = rel_prop, color = primaryDiagnosis)) +
  geom_point(alpha = 0.5, size = 0.5) + 
  geom_smooth(se = F, method = 'lm', fullrange = T) +
  ylab('Relative Cell Prop.') + 
  xlab('Age') + 
  theme(
    axis.text.x = element_text(size = 15),
    axis.text.y = element_text(size = 15),
    axis.title.x = element_text(vjust = -2, size = 15),
    axis.title.y = element_text(size = 15),
    legend.position = c(0.9, 0.9),
    legend.title = element_blank(),
    legend.text = element_text(size = 12),  
    strip.text.x = element_text(size = 14)) +
  facet_grid(~cell_type, drop = T, scale = "free_x", space = "free") +
  guides(fill = guide_legend(nrow = 1)) +
  scale_y_continuous(limits = c(-3, 3)) + 
  scale_color_manual(values = c("dodgerblue2", "firebrick2"), labels=c('CON', 'SCZ'))

Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2 3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

#Figure 2d: Mega-analysis results when binned as +/- 70 years
significant_mega_age_class = subset(mega_lms_by_age_class, padj <= 0.1 & term == "SCZ" & cell_type %in% c("SST", "PVALB", "VIP"))
significant_mega_age_class$asterisks = ifelse(significant_mega_age_class$padj <= 0.01, "***", ifelse(significant_mega_age_class$padj <= 0.05, "**", "*"))
significant_mega_age_class$vjust = ifelse(significant_mega_age_class$estimate >= 0, -3, 3)

figure_2d = mega_lms_by_age_class %>% 
  filter(term %in% 'SCZ' & cell_type %in% c("SST", "PVALB", "VIP")) %>% 
  mutate(cell_type = fct_reorder(cell_type, estimate)) %>% 
  ggplot(aes(x = age_class, y = estimate, fill = age_class)) + 
  geom_hline(yintercept = 0) + 
  geom_bar(stat = "identity", position = "dodge") + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error)) + 
  ylab('Beta Coefficient') + 
  xlab('\n \n Age Range') + 
  scale_fill_manual(values = c("mediumpurple3", "olivedrab3"), labels=c('Above 70', 'Below 70')) +
  geom_text(data = significant_mega_age_class,
            aes(age_class, estimate, label = asterisks), 
            vjust = significant_mega_age_class$vjust) +
  theme(
    axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1, size = 13), 
    axis.text.y = element_text(size = 15),
    axis.title.x = element_text(vjust = 8, size = 15),
    axis.title.y = element_text(size = 15),
    strip.text.x = element_text(size = 14)) +
  facet_wrap(~cell_type, drop = T, scale = "free_x") +
  guides(fill = FALSE) +
  scale_y_continuous(limits = c(-0.55, 0.55))

figure_2_middle = (figure_2b +figure_2d) + plot_layout(widths = c(2,1))

figure_2 =  figure_2c / figure_2_middle / figure_2a
figure_2 = figure_2 + plot_annotation(tag_levels = 'A') + plot_layout(heights = c(1, 1, 1)) & theme(plot.tag = element_text(size = 25))
figure_2

`geom_smooth()` using formula = 'y ~ x'

Warning: Removed 53 rows containing non-finite outside the scale range (`stat_smooth()`).
Warning: Removed 53 rows containing missing values or values outside the scale range (`geom_point()`).

#Save to plots folder
output_filename <- "~/cell_prop_psychiatry/plots/figure_2.jpg"

# Use ggsave to save the plot as a JPG image
ggsave(output_filename, figure_2, width = 13, height = 15, dpi = 500)

`geom_smooth()` using formula = 'y ~ x'

Warning: Removed 53 rows containing non-finite outside the scale range (`stat_smooth()`).
Warning: Removed 53 rows containing missing values or values outside the scale range (`geom_point()`).

Figure 3: Validating trends and interactions in ssRNAseq

# Figure 3a: snRNAseq cell proportions binned by age
significant_sn_mega_by_age <- subset(sn_mega_lms_by_age, padj <= 0.1 & term == "SCZ" & cell_type %in% c("SST", "PVALB", "VIP"))
significant_sn_mega_by_age$asterisks <- ifelse(significant_sn_mega_by_age$padj <= 0.01, "***",
                                         ifelse(significant_sn_mega_by_age$padj <= 0.05, "**", "*"))
significant_sn_mega_by_age$vjust <- ifelse(significant_sn_mega_by_age$estimate >= 0, -5, 5)

figure_3a <- sn_mega_lms_by_age %>%
  filter(term %in% 'SCZ' & cell_type %in% c("SST", "PVALB", "VIP")) %>%
  mutate(cell_type = fct_reorder(cell_type, estimate)) %>%
  ggplot(aes(x = age_group, y = estimate, fill = age_group)) +
  geom_hline(yintercept = 0) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error)) +
  ylab('Beta Coefficient') +
  xlab('Age Range (Above/Below 70)') +
  scale_fill_manual(values = c("mediumpurple3", "olivedrab3"), labels = c('Above 70', 'Below 70')) +
  geom_text(data = significant_sn_mega_by_age,
            aes(age_group, estimate, label = asterisks),
            vjust = significant_sn_mega_by_age$vjust) +
  theme(
    axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1, size = 13),
    axis.text.y = element_text(size = 15),
    axis.title.x = element_text(vjust = 1, size = 15),
    axis.title.y = element_text(size = 15),
    strip.text.x = element_text(size = 14),
    legend.position = "none") +
  scale_y_continuous(limits = c(-1, 1)) +
  facet_wrap(~cell_type, drop = TRUE, scale = "free")

# Figure 3b: Relative cell proportions from snCTPs
figure_3b <- snCTPs_long %>% filter(total_cells >= 500) %>%
  filter(Phenotype %in% c('SZ', "CON"), cell_type %in% c("SST", "PVALB", "VIP")) %>%
  ggplot(aes(x = Age, y = rel_prop, color = Phenotype)) +
  geom_point(alpha = 0.8, size = 1.3) + 
  geom_smooth(se = F, method = 'lm') +
  scale_colour_hue(labels = c("Control", "SCZ")) +
  ylab('Relative Cell Prop.') + 
  xlab('Age') + 
  theme(
    axis.text.x = element_text(size = 13),
    axis.text.y = element_text(size = 13),
    axis.title.x = element_text(, size = 16),
    axis.title.y = element_text(size = 16),
    legend.position = c(0.94, 0.94),
    legend.title = element_blank(),
    legend.text = element_text(size = 12), 
    strip.text.x = element_text(size = 15),
    strip.text.y = element_text(size = 15)
  ) +
  facet_rep_grid(rows = vars(Institution),
    cols = vars(cell_type),
    drop = T,
    scale = "free",
    switch = "y",
    space = "free_x",
    repeat.tick.labels = T
  ) +
  guides(fill = guide_legend(nrow = 1)) + 
  scale_color_manual(values = c("dodgerblue2", "firebrick2"), labels=c('CON', 'SCZ'))

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

# Figure 3c: Interaction terms for snRNAseq results
significant_sn_mega <- subset(sn_mega_lms, padj <= 0.1 & term == "Interaction" & cell_type %in% c("Astrocyte", "IT", "L5.6.NP", "LAMP5", "Microglia", "Oligodendrocyte", "OPC","Pericyte", "PVALB", "SST", "VIP" ))
significant_sn_mega$asterisks <- ifelse(significant_sn_mega$padj <= 0.01, "***",
                                         ifelse(significant_sn_mega$padj <= 0.05, "**", "*"))
significant_sn_mega$vjust <- ifelse(significant_sn_mega$estimate >= 0, -7, 7)

sn_mega_lms$cell_type = sn_mega_lms$cell_type %>% factor(levels = c("PVALB", "SST", "VIP", "LAMP5", "PAX6", "IT", "L4.IT", "L5.ET", "L5.6.IT.Car3", "L5.6.NP", "L6.CT", "L6b", "Astrocyte", "Endothelial", "Microglia", "Oligodendrocyte", "OPC", "Pericyte", "VLMC"))
figure_3c = sn_mega_lms %>% filter(cell_type %in% c("Astrocyte", "IT", "L5.6.NP", "LAMP5", "Microglia", "Pericyte", "PVALB", "SST", "VIP" )) %>%
  filter(term %in% 'Interaction') %>% 
  ggplot(aes(x = cell_type, y = estimate)) + 
  geom_hline(yintercept = 0) + 
  geom_bar(stat = "identity", position = "dodge") +
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error)) + 
    facet_grid(~cell_class, drop = T, scale = "free_x", space = "free") +
  ylab('Interaction Coefficient') + 
  xlab('Cell Type') + 
  geom_text(data = significant_sn_mega,
            aes(cell_type, estimate, label = asterisks), 
            vjust = significant_sn_mega$vjust) +
  theme(
    axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1, size = 12), 
    axis.text.y = element_text(size = 15),
    axis.title.x = element_text(vjust = -2, size = 15),
    axis.title.y = element_text(size = 15),
    strip.text.x = element_text(size = 15)) +
  guides(fill = FALSE) +
  scale_y_continuous(limits = c(-0.45, 0.65))


# Figure 3: Combined figure
figure_3_top = figure_3b
figure_3_bottom  = figure_3c + figure_3a + plot_layout(widths = c(2,1)) 

figure_3 <- figure_3_top/figure_3_bottom

figure_3 = figure_3 + plot_annotation(tag_levels = 'A') + plot_layout(heights = c(2,1)) & theme(plot.tag = element_text(size = 25))

figure_3

`geom_smooth()` using formula = 'y ~ x'

#Save to plots folder
output_filename <- "~/cell_prop_psychiatry/plots/figure_3.jpg"

# Use ggsave to save the plot as a JPG image
ggsave(output_filename, figure_3, width = 13, height = 16, dpi = 500)

`geom_smooth()` using formula = 'y ~ x'

Figure 5: AgeOnset vs. MGP

#Duration of Illness
mgp_estimations_long$durationIllness <- mgp_estimations_long$ageDeath - mgp_estimations_long$ageOnset

CELL_TYPES = c("SST", "PVALB", "VIP")
DATASETS = c("GVEX", "LIBD_szControl")
ageOnset_residuals_data = matrix(ncol = (ncol(mgp_estimations_long)+1), nrow = 0) %>% data.frame()
colnames(ageOnset_residuals_data) = c(colnames(mgp_estimations_long), "resid")

for(DATASET in DATASETS) {
  for(CELL_TYPE in CELL_TYPES) {
    #Create res_model without ageOnset variable
    res_model_data = mgp_estimations_long %>% filter(primaryDiagnosis %in% c('Schizophrenia'), cell_type == CELL_TYPE , dataset == DATASET, is.na(ageOnset) == F)
    #Model without ageOnset
    res_model = lm(scale(rel_prop) ~ scale(ageDeath) + scale(RIN) + scale(PMI) + reportedGender, data = res_model_data)
    #Calculate residuals and add column to res_model_data
    res_model_data = res_model_data %>% add_residuals(var = "resid", model = res_model)
    ageOnset_residuals_data = ageOnset_residuals_data %>% rbind(res_model_data)
  }
}

figure_5 <- ageOnset_residuals_data %>%
  ggplot(aes(x = ageOnset, y = resid, color = dataset)) +
  geom_point(alpha = 0.5, size = 0.8) + 
  geom_smooth(se = FALSE, method = 'lm', fullrange = TRUE) +
  ylab('rCTP residuals (AU)') + 
  xlab('Age of SCZ onset') + 
  theme(
    axis.text.x = element_text(size = 16),
    axis.text.y = element_text(size = 16),
    axis.title.x = element_text(vjust = -2, size = 16, margin = margin(t = 8)),
    axis.title.y = element_text(size = 16),
    strip.text = element_text(size = 16),
    legend.text = element_text(size = 16),
    plot.margin = margin(10, 10, 20, 10)) +
  facet_grid(cols = vars(cell_type), rows = vars(dataset)) +
  scale_y_continuous(limits = c(-1.5, 1.5)) + 
  scale_color_manual(values = c("dodgerblue2", "firebrick2")) +
  guides(color = guide_legend(title = NULL)) +
  stat_cor(aes(label = paste(..r.label.., ..p.label.., sep = "~")), color = "black", geom = "label")


figure_5

Warning: The dot-dot notation (`..r.label..`) was deprecated in ggplot2 3.4.0.
ℹ Please use `after_stat(r.label)` instead.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

`geom_smooth()` using formula = 'y ~ x'

Warning: Removed 34 rows containing non-finite outside the scale range (`stat_smooth()`).
Warning: Removed 34 rows containing non-finite outside the scale range (`stat_cor()`).
Warning: Removed 34 rows containing missing values or values outside the scale range (`geom_point()`).

#Save to plots folder
output_filename <- "~/cell_prop_psychiatry/plots/figure_5.jpg"

# Use ggsave to save the plot as a JPG image
ggsave(output_filename, figure_5, width = 9, height = 7, dpi = 500)

`geom_smooth()` using formula = 'y ~ x'

Warning: Removed 34 rows containing non-finite outside the scale range (`stat_smooth()`).
Warning: Removed 34 rows containing non-finite outside the scale range (`stat_cor()`).
Warning: Removed 34 rows containing missing values or values outside the scale range (`geom_point()`).